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We study the equilibrium time correlations for the conserved fields of classical 
anharmonic chains and argue that their dynamic correlator can be predicted on the 
basis of nonlinear fluctuating hydrodynamics. In fact our scheme is more general and 
would cover also other one-dimensional hamiltonian systems, for example classical 
and quantum fluids. The only input parameters required are the average equilibrium 
currents and the static susceptibilities of the conserved fields. In our context fluc- 
tuating hydrodynamics is a nonlinear system of conservation laws with noise. For a 
single mode it is equivalent to the noisy Burgers equation, for which explicit solutions 
are available. Our focus is the case of several modes. No exact solutions have been 
found so far and we rely on a one-loop approximation. The resulting mode-coupling 
equation has a quadratic memory kernel and describes the time evolving correlator 
matrix of all locally conserved fields. Long time asymptotics is computed analyti- 
cally and finite time properties are obtained through a numerical simulation of the 
mode-coupling equations, which provide predictions still to be compared with data 
from molecular dynamics. 
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To predict the dynamic correlator of anharmonic chains is still a theoretical challenge. 
In higher dimensions fluctuating hydrodynamics serves as a convenient starting point [1-3] . 
But, as recognized already in the 1970ies [4, 5], in one dimension, while the static cor- 
relations are of short range, the dynamic current-current correlations generically have an 
anomalously slow decay. In particular the transport coefficients, required as an input for 
fluctuating hydrodynamics, are no longer properly defined. There have been huge efforts, 
both through theoretical investigations and numerical simulations, to more precisely char- 
acterize this "anomalous" behavior (for a partial list on FPU chains only, see [6-11]). Here 
we argue that, in one dimension, linear fluctuating hydrodynamics has to be extended to a 
nonlinear version, which will be outlined below. There have been related attempts before 
[12, 13]. Our main advance is to treat the full system of coupled conserved modes and to 
run time- dependent numerical simulations of the respective mode-coupling equations. In 
the simulation we use the exact, microscopically computed parameters for the particular 
anharmonic chain under consideration and thereby provide time-resolved predictions which 
can be compared directly with the corresponding molecular dynamics. 

As a start-up, let us recall the case of a single conserved field, say u(x,t), space x G E, 
time t, which satisfies the conservation law 



with given current function j(u). We want to study the fluctuations relative to a uniform 
background u, i.e. u(x, t) = u+u(x, t), hence expand (1) to second order and add dissipation 
and noise, 



where £,(x,t) is normalized space-time white noise. We consider the mean zero, space- 
time stationary process u(x,t) governed by (2). Then the spatial statistics at fixed time 
t is white noise, (u(x,t)u(x' , £)) = x${ x ~ x ')i X — o~/2D, which reflects that the static 
correlations of an underlying microscopic model decay exponentially fast. Of particular 
interest is the correlator S(x,t) = (u(x,t)u(0,0)), S(x,0) = X$( x )- Its large scale behavior 
will be dominated by the nonlinearity, but dissipation and noise is required to maintain the 
proper steady state. (2) is the noisy Burgers equation, equivalently the spatial derivative 
of the one-dimensional KPZ equation [14]. There is an exact computation of S(x,t) using 




(1) 



8 t u + d x (f(u)u + \]"{u)u 2 - Dd x u + V^O = , 



(2) 
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replica [15]. In particular one knows the universal long time limit, 

5(x,t) = X (AB|t|)^ 2/3 /KPz((A B |t|)- 2/3 (a:-j , (u)t)), (3) 

valid for large x,t, with A| = 2xj"{u) 2 . Identical scaling properties have been derived also 
for stochastic lattice gases [16, 17]. The universal scaling function /kpz can be written in 
terms of a Fredholm determinant and has been computed with great precision [18]. Because 
of the nonlinearity the spreading is faster than diffusive. On the other hand D, a appear 
in (3) only through the static covariance x- This quantity will have to come from the 
microscopic model, while dissipation and noise by themselves are phenomenological, not 
directly identifiable with microscopic properties. (3) and similar type of predictions based 
on (2) have been confirmed experimentally using a thin film of turbulent liquid crystal 
[19] and in Monte Carlo simulations of Eden cluster growth [20]. These results provide an 
indirect confirmation that higher order nonlinearities, like u 3 , will not change the scaling 
asymptotics. 

Anharmonic chains have three conservation laws, compression (or elongation), momen- 
tum, and energy. Thus we have to extend the previous considerations to several modes (or 
components). We use a as mode index. Then (1) generalizes to 

d t u a + d x j a (u) = , a = 1, . . . , n , (4) 

u — (ui, ... , u n ). Expanding as u a = u a + u a , the coefficients of the linearized equation are 

A afS (u) = d Uf J a (u) (5) 

and coefficients of the quadratic part are given by the Hessians 

#a 7 (if) = d Up d u J a (u) . (6) 

Since the background u is fixed, it will be suppressed in our notation. The equal time 
equilibrium correlations decay rapidly. Hence u a (x,t) at fixed t is modeled as white noise 
with covariance (u a (x,t)up(x',t)) = C a p5(x — x'). C Q p = Cp a and C > as a matrix. As 
only microscopic input, nonlinear hydrodynamics requires the currents j a , more precisely 
A,H a , and the static correlator C. As can be verified directly, for anharmonic chains it 
holds 

AC = CA l (7) 
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with 't' denoting transpose. Hence A has a system of right and left eigenvectors, Aip a = 
c a ip a , A t i/) a = c a ip a , and it is convenient to switch to normal modes. We form the matrices 

R = KVU){o=i,„.,n} ,R l = (l^o)) {a=i,...,n}- U P to a g lobal si g n tne eigenvectors are unique 
when normalized such that RR~ l = 1 , .RCi?* = 1. Then RAR~ X = diag(ci, . . . , c n ) and 
the normal modes, (p a , are defined through <ft = Au. While not strictly necessary, the 
linear transformation to normal modes is of advantage, since normal modes have a definite 
propagation velocity and turn out to be essentially uncorrelated. 

We now expand (4) to second order in u, transform to normal modes, and add dissipation 
and noise with the result 

n 

d t K + d x (c a <p Q + G%M-r - 9 x {D<f>) a + Q=0, (8) 

where = | Y^ a ' R aa' p\H a ' ^•y) ■ The diffusion matrix D is positive definite. £ a (x,t) is 
space-time white noise with strength 

(Ux, t%{x', t')) = 2D afi 5(x - x')5(t - t') . (9) 

As before, we consider the mean zero, stationary process 4>(x,t) governed by (8). In the 
linear case, G = 0, 4>(x, t) is a Gaussian process, which for fixed t has white noise statistics 
with independent components of unit strength as required by RCR 1 — 1. It can be checked 
that, for a suitable spatial discretization of (8), this measure remains invariant under the 
full nonlinear evolution, provided the couplings are symmetric in all indices [21]. In 
general, we assume that the steady state of (8) has short range correlations. 
On the basis of (8) one would like to compute the dynamic correlator 

S aP (x,t) = (0 a (M)<MO,O)>, (10) 

which would serve then as an approximation to the true dynamic correlator of the micro- 
scopic model. Because of the nonlinearity, this is a difficult program and we rely here on 
a mode-coupling equation resulting from the one-loop expansion of (8). The precise details 
of the computation can be found elsewhere [21]. In essence, one writes down the evolution 
equation for S a p{x,i), once iterated in time, which yields a time integral over the Gaussian 
propagator and an expectation of the form (0(s)0(s)0(O)0(O)) for the intermediate time s, 
< s < t. The Gaussian propagator is replaced by the full propagator and the 4-point 
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correlation is replaced by its Gaussian pairings, where the equal time pairing vanishes be- 
cause of the conservation law and the remaining two pairings yield the same contribution. 
The result is a cubic memory equation. Anticipating already the application to anharmonic 
chains with integer particle index, numerically we use a unit spatial grid. Adopting the 
standard conventions for discrete Fourier transforms, in Fourier space the mode-coupling 
equation reads 

d t S a p(k,t) = -ic a sin(2nk)S a p(k,t) 

n „i 

- 2(1 - cos(27rA;)) V (-D aa ,S a ,p{k,t)+ / dsM aa ,(k,s)S a/ p(k,t-s)), (11) 

a'=l Jo 
k E [— |, |], with the quadratic memory kernel 

M aa ,(k,t) = 2 G %r G P'f / <k$wfr ~ Q,t)S 77 ,(q,t) , (12) 

^,/3',7 )7 '=l J -1/2 

which has to be solved with the initial condition 

S> a p{k,O) = 6 a 0. (13) 

(11) is a lowest order approximation. However, for n = 1, first written down in [22], the 
exact solution (2) can be compared with the numerical solution of (11), (12). In Fig. 1 we 
display a time sequence for a single mode, n = 1, with G\ x = |. For t > 32 the scaled 
solution remains stationary. We compare with the universal KPZ scaling function, thereby 
reproducing Fig. 3 of [23]. The asymptotic scaling function differs from /kpz by a few 
percent only. We expect that such a precision extends to several modes. 

Numerically one observes that, while off-diagonal elements of S a p develop immediately, 
they decay fairly rapidly. Thus we may invoke the diagonal approximation, S a p(x,t) ~ 
5 a pf a {x,t) in position space, and the memory kernel becomes 

n 

M aa (x,t) = 2 J2(G a ^) 2 fp(x,t)f 7 (x,t). (14) 

Let us assume that all velocities {c a } are distinct. Then, ignoring nonlinearities, the mode 
f a is sharply peaked at c a t and the modes separate from each other linearly in t. Hence the 
product fpfry decays rapidly, unless j3 = 7, and the off-diagonal terms in the sum (14) can 
be discarded. In general, two distinct scenarios may occur, (i) G^ a 7^ 0. Then G^afafa 
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FIG. 1. (Color online) Time sequence of the numerical mode-coupling solution, n = 1, = ^, 
rescaled by (\/2t) 2 / 3 . The curves are piecewise-linear due to the underlying spatial lattice. The 
exact KPZ scaling function is shown in dashed red. 

travels with the same velocity as f a and dominates the interaction with the other modes. 
The dynamics of the mode f a reduces to the case n = 1 with scaling (2) and shift c a t. This 
is the KPZ dominated behavior, (ii) G° a = 0. The leading term vanishes and long time 
behavior is dominated by ^ 0. If also ^ 0, then the coupling of mode a to the 
pair of modes (3 results in a scaling function f a (k,t) = exp[— |/c| 5//3 |Aht|] with computable 
coefficient Ah, see [21]. 

There are more cases to be considered, but equipped with this nonlinear extension of 
fluctuating hydrodynamics it is more interesting to understand how one obtains a prediction 
for the structure factor of anharmonic chains. We first have to figure out the conserved 
fields and their macroscopic Euler equations. The chain consists of N particles, position qj, 
momentum pj, j = 1, . . . , N, unit mass, and is governed by the hamiltonian 

N 

i=i 

where periodic boundary conditions, qN+i — q± + L, are imposed. A prototypical potential 
is the FPU choice V(y) = |y 2 |ay 3 + \y A . The locally conserved microscopic fields are 
compression rj = g J+1 — qj, momentum pj, and energy ej = \p 2 + Vijj). In a microcanonical 
simulation one fixes the length (compression) per particle, £, as L = N£, the momentum per 
particle, u, as ^=iPj = Nu, an d the energy per particle, e, as H = Nt. Computationally, 
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it is convenient to switch to the grand canonical pressure ensemble. Then £ is paired with 
the pressure p and the internal energy e with the inverse temperature /3. In the grand 
canonical ensemble, {rj,pj} become independent random variables. The distribution of pj is 
a Maxwellian shifted by u, and the distribution of rj is given by Z~ x exp[—(3(V(y) + py)\ = 
(•) p ,p with partition function Z = J dyexp[— (3(V(y) + py)}. Clearly the pressure is just 
the average force acting on a specified particle. The microcanonical and grand canonical 
parameters are related through 

t = <v)pj>, e=^ + (V(y))p,P- (16) 
The hydrodynamic currents are 

jt = -u , j u = p(i, e - ^u 2 ) , j t = up(£, c - \u 2 ) , (17) 

which inserted in (4) result in the Euler hydrodynamics of the anharmonic chain. Note that 
variable x corresponds to the continuum approximation of the particle index. 

Without loss of generality we expand at u = 0. All coupling coefficients appearing in 
(12) can be computed in terms of at most third order cumulants involving y,V(y). These 
integrals and the somewhat unwieldy substitutions are easily performed using Mathematica. 
There are three modes: the heat mode, a = 0, with velocity Co = and two sound modes, 
a = ±1, with velocity c a = ac, a = ±1, where c is the sound speed, 

c 2 = — dtp + pd e p . (18) 

Since momentum is conserved, H a o = 0. Most importantly G® = 0, always. Hence the 
heat mode is predicted to have non-KPZ scaling. A particular and much studied case is 
an even potential, V(y) = V(—y), at p = 0. Then many cumulants vanish. In fact the 
only non-zero coefficients are Gq ± and G°^ + ,G°__. On this basis one expects the sound 
modes to be diffusive and the feed-back on the heat mode to lead to the scaling function 
/ (fc,t)=exp[-|fc| 3 / 2 A h |t|]. 

To go beyond the large t asymptotics we turn to the simulation of the mode-coupling 
equations, which is based on the Fourier space representation (11), in such a way that the 
values of the memory kernel M aa >(k,s) for s < t can be stored and re-used. The time 
and momentum variables are discretized by a uniform grid. The parameters of the chain 
are picked somewhat arbitrarily but related to the parameters used in molecular dynamics 
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FIG. 2. (Color online) Time sequence of normal mode correlations for the FPU chain with a = 2, 
p = 1, (3 = 0, and u = (a), (b), (c). Magnification of the central heat mode peak (d) and the 
right sound mode peak (e) in suitably rescaled coordinates. 

[9, 10]. We stress that simulations can be performed for any choice of the potential and 
thermodynamic parameters at minimal numerical efforts. Fig. 2 is obtained for the FPU 
potential with a = 2 and u = 0, p — 1, j5 = 2, resulting in c = 1.455. V{y) + py has a single 
critical point at y = —1.755. The grey vertical lines at ±ct indicate the predicted position 
of the sound mode peaks. In the time sequence we display the superimposed normal mode 
correlations (area 1 under each curve). More details are provided in the blow-up. For the 
heat mode one observes oscillations which move away from the center and eventually die out. 
The tail of the heat mode is cut off at the location of the sound mode, but in-between the two 
sound mode peaks the theoretically predicted curve is approached. At the longest available 
time the sound modes are still asymmetric and have not yet reached their asymptotic shape. 

Based on these and further simulations of the mode-coupling equations for anharmonic 
chains, the following qualitative picture for the motion of the peaks in index number space 
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seems to emerge. The sound modes "rapidly" decay to a shape function which is centered at 
act and varies on the scale t 2//3 . The shape function itself is still slowly varying. The coupling 
G° aa determines the scaling of the heat mode. Since only the integral over the square of the 
shape function is involved, the heat mode rapidly achieves its asymptotic form in the range 
{\x\ < ct} and with a still slowly varying non-universal constant. The slow motion of the 
sound modes is governed by Gq and G a _ a _ a . Assuming already the overall scaling picture, 
the size of these corrections is estimated to be of the order t -1 / 15 , resp. t -1//9 , relative to 
the leading term which indicates that /kpz is approached rather slowly. Of course, only a 
guideline can be presented. For the precise dynamics all couplings have to be used. 
Conclusions. We developed a nonlinear extension of fluctuating hydrodynamics applicable 
to one-dimensional systems, in principle including classical fluids, quantum fluids [24-26], 
and quantum spin chains. Already at the level of the one-loop approximation it is crucial 
to maintain the couplings between all conserved modes. Numerical solutions of the mode- 
coupling equations provide a realistic picture of the correlation dynamics. On the limited 
time scale available, the numerical solutions are in agreement with the analytical computa- 
tions. But the theory still has to pass the key test, namely a quantitative comparison with 
molecular dynamics. In this test one has to keep in mind that nonlinear fluctuating hydro- 
dynamics is itself an approximation at best valid only for long times and wavelengths. In 
addition mode-coupling is an approximation to nonlinear fluctuating hydrodynamics, albeit 
likely to be fairly accurate, as demonstrated in Fig. 1 
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